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NUMERICAL SOLUTION OF THE UNSTEADY NAVIER -STOKES EQUATIONS 
AND APPLICATION TO FLOW IN A RECTANGULAR CAVITY 
WITH A MOVING WALL 
by Leo F. Donovan 
Lewis Research Center 

SUMMARY 

A computer program to solve the unsteady, two-dimensional, incompressible 
Navier-Stokes equations was written in FORTRAN IV. The numerical method makes use 
of an iterative solution of a Poisson's equation for pressure followed by an explicit cal- 
culation of velocities. The computer program is included as an appendix. 

Unsteady flow in a two-dimensional, rectangular cavity with the upper wall moving 
at constant velocity is investigated using the computer program. The calculations start 
with the fluid at rest in the cavity and continue until no further changes in velocity occur. 
Results for cavities with aspect ratios of 1/2, 1, and 2 are presented for a Reynolds 
number of 100. For a square cavity, results are also given for several Reynolds num- 
bers between 100 and 500. Velocities calculated from the unsteady Navier-Stokes equa- 
tions at large times are compared where possible to velocities calculated from the steady 
Navier-Stokes equations and to the results of steady experiments; good agreement is 
shown. 

A technique for conducting a numerical flow visualization experiment in conjunction 
with the solution of the Navier-Stokes equations is described. The results of the experi- 
ment are recorded on film which can be shown as a motion picture. Selected frames 
from the motion picture made during the investigation of cavity flow are reproduced in 
this report to illustrate this method of data presentation. 


INTRODUCTION 

The availability of large, high-speed computers allows us to attack some formidable 
and interesting problems. One such area of research is in the solution of nonlinear par- 
tial differential equations describing physical phenomena. Analytic solutions can only be 



obtained for certain special cases. In addition, the nonlinearity of the equations makes 
them much more difficult to solve numerically than linear equations. However, some 
progress has been made. The book by Ames (ref. 1) summarizes much of this work. 

An additional difficulty with numerical solutions of partial differential equations is 
that they provide a mass of detailed information that is very hard to assimulate and 
understand. In order to overcome this problem Fromm and Harlow (ref. 2) have devised 
a visual display technique for presenting the results of their fluid dynamics calculations. 
This technique is analogous to a flow visualization experiment in the laboratory, in which 
a tracer is introduced into a fluid to make the flow visible. 

Fluid motion is governed by the continuity equation and Navier-Stokes equations, 
expressing conservation of mass and momentum. Almost all the numerical solutions of 
these equations have been for two-dimensional flows. Some investigators have trans- 
formed the equations to stream function and vorticity coordinates. Pearson (ref. 3), for 
example, has treated rotating disks in this manner. Other investigators have chosen to 
retain velocity and position coordinates. Work by Harlow and coworkers (ref. 4) on free 
surface flows falls in this category. Since Harlow's method has the advantages of com- 
bining a successful numerical technique with visual display we have used it without the 
free-surface feature in our studies. 

The results of a numerical investigation of incompressible flow in a square cavity at 
a Reynolds number of 100 are presented in reference 5. The present report describes 
an extension of the technique to rectangular cavities and higher Reynolds numbers. The 
differential equations describing unsteady, incompressible flow are discussed first. The 
corresponding difference equations are then derived and the numerical method used to 
solve them is presented. The method is then used to calculate the startup of flow in a 
rectangular cavity with a moving wall. 

An advantage of the technique is that unsteady results are obtained. However, since 
no unsteady experimental results or prior calculation for cavity flow are available, only 
comparisons at steady conditions are possible. Results for rectangular cavities with 
aspect ratios of 1/2 and 2 are presented for a Reynolds number of 100. For flow in a 
square cavity, results are given for Reynolds numbers from 100 to 500. 


ANALYSIS 

Differential Equations 

Conservation of mass and momentum are expressed by the continuity equation and 
Navier-Stokes equations, respectively. For constant density and viscosity the two- 
dimensional forms of these equations for a Newtonian fluid are (ref. 6) 
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these equations u and v are the velocity components in the x and y directions, 

is the pressure, and g and g,, are the body forces in the x and y directions. 

x y ■ 

The overbars are used to denote dimensional quantities. 

The equations can be made dimensionless with a reference velocity W, a reference 
length L, and the fluid viscosity v by the following substitutions: 
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The viscosity v is thus the reciprocal of the Reynolds number. 
The dimensionless forms of equations (1) to (3) then become 
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The continuity equation can be used to write equations (6) and (7) in forms such that 
the resulting finite difference equations will rigorously conserve momentum (see ref. 4). 
Thus 
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If the x- and y-momentum equations are differentiated with respect to x and y, 
respectively, and the results are added and rearranged, the following Poisson's equation 
for pressure is obtained 
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The first term to the right of the equality sign in this equation is the time derivative of 
the left side of the continuity equation and, as such, should be zero. However, since the 
continuity equation will not be satisfied exactly, this term is retained as a correction. 


Finite Difference Equations 

Computational mesh . - The finite difference mesh is shown in figure 1. The posi- 
tions of the variables are chosen so that vertical walls pass through positions where the 
u- component of velocity is defined, horizontal walls pass through positions where the 
v-component of velocity is defined, and pressures are cell-centered. These positions 
have been chosen to facilitate application of the boundary conditions. 

If it is necessary to evaluate one of the variables at a position where it is not de- 
fined, an average is used. For example, 
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Figure 1. - Computational mesh. 


For conciseness, however, terms such as u. , will be retained in the equations through- 
out the remainder of this section. 

Representation of derivatives . - When converting the differential equations to finite 

difference form, centered differences are used to represent derivatives. For example, 
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When the derivative of a product of undefined quantities is formed, the product is differ- 
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entiated and then averages are formed. For example, 9 (uv)/9x 9y at the point (i,j) is 
replaced by 
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A typical average is then 

(uv) i+l/2, j+1/2 = \ (u i+l, j+1/2 + u i, j+l/2 )(v i+l/2, j+1 + v i+l/2, (12) 


Difference equations . - If the abbreviation 

d i,j = ^ (Ui ,j+i/2 " U i,j-l/2 } + ^ (v i+l/2,j " v i-l/2,j ) 
is introduced, the finite difference form of the continuity equation can be written as 
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Formally, the time derivative of the continuity equation is then 
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where the superscripts n + 1 and n refer to the advanced time and the current time, 
respectively. Hereinafter, the lack of a superscipt will indicate the current time. 

The finite difference Poisson's equation for pressure can be obtained from equa- 
tion (10) 
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The final term arises because equation (16) will be solved iteratively and therefore the 

continuity equation will not be satisfied exactly. The discrepancy d” . from the cur- 

J 

rent time is used as a correction during the calculations at the advanced time. Since it 
is desired that the continuity equation be satisfied as closely as possible at the advanced 
time, the d^ + .^ are set to zero. If the correction terms were not included in equa- 
tion (16), the pressure iteration would have to be carried further; this would require 
more computer time than the technique used here. Hirt and Harlow (ref. 7) discuss the 
use of a correction term in more detail. 

The finite difference forms of equations (2) and (3) for the velocities are 
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These two equations are solved at each time step after the pressure has been determined. 

Initial and Boundary Conditions 

The specific problem to be solved will determine the initial and boundary conditions. 
At a no-slip wall the tangential component of velocity is equal to the wall velocity; if the 
wall is impermeable, the normal component of velocity is zero. At slip walls (i. e. , 
planes of symmetry) the normal component of velocity is zero and the normal gradient of 
the tangential component is zero. If fluid is entering and leaving the computing region, 
sufficient information must be provided at the inlet and outlet to make the problem deter- 
minant. If the fluid starts from rest, the initial conditions are simply that the velocities 
are everywhere zero and the pressure is uniform. 

For cavity flow, which is discussed later in this report, the appropriate boundary 
conditions are no-slip, impermeable walls. One of the walls is moving at constant ve- 
locity and the other three walls are stationary. 

No-slip wall . - Fictitious tangential velocities outside the cavity are defined because 
they are needed when equations (16) to (19) are solved for the row of cells just inside the 
cavity. These fictitious tangential velocities are defined so that the interpolated values 
of velocity at the wall satisfy the no-slip boundary condition. The fictitious velocities 
are calculated at each time step, after the velocities inside the cavity have been evalu- 
ated. 

For a vertical wall at j - (1/2), figure 2(a), the fictitious velocity v.^jyg j_i 
would be evaluated so that 
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where the constant is the wall velocity, which is zero for a stationary wall or unity for 
the moving wall. 

Similarly, for a horizontal wall at i - (1/2), figure 2(b), the fictitious velocity 

u. . . 1 /0 would be evaluated so that 
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where the constant is again the wall velocity. 
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Figure 2. - Boundary conditions. 


I mpermeable wall . - The normal velocity at an impermeable wall can be written 
directly. For a vertical wall at j - (1/2), figure 2(a), 

Uj ~ ® * (22) 

For a horizontal wall at i - (1/2), figure 2(b), 

v i-l/2, j = 0 a11 i (23) 

Pressure . - Fictitious pressures outside the cavity are also defined because they 
are needed in the solution of equation (16) for the row of cells just inside the cavity. The 
Navier-Stokes equations can be evaluated at the walls to determine the fictitious pres- 
sures since the normal component of velocity is always zero at an impermeable wall. 

For a left stationary wall at j - (1/2), figure 2(a), the fictitious pressure cp ■ ■ •. can 
be calculated from equation (18) as 
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For a stationary bottom wall at i - (1/2), figure 2(b), the fictitious pressure <p. . can 
be calculated from equation (19) as 


%] = ^i-lj (u i-l, j+1/2 ’ u i-l,j-l/2 } 


fix 


(25) 


The velocities are known, either from the initial condition or the previous time cycle. 


Numerical Flow Visualization Experiment 

The traditional way of presenting the results of experimental or numerical fluid 
mechanical investigations is not the best way of helping one form a coherent, overall 
picture of a complicated flow situation, especially if the flow is unsteady. Flow visual- 
ization experiments, in which a tracer is introduced into the fluid to make its movement 
visible, have been designed to overcome this difficulty in the laboratory. A numerical 
analog of a laboratory flow visualization experiment offers the same advantages for nu- 
merical fluid dynamics studies. The technique devised by Fromm and Harlow (ref. 2) 
employs special marked particles that move with the fluid as the tracer. A microfilm 
recorder is used to photograph the marked particles displayed on a cathode ray tube. 

The sequence of photographs, when viewed as a motion picture, shows the behavior of 
the fluid clearly. The marked particles serve only to make the flow visible and do not 
enter into the solution of the Navier-Stokes equations. 

The flow visualization experiment is conducted as follows: At the start of the calcu- 
lation an initial uniform distribution of either one or four particles per cell is imagined 
to be distributed throughout the fluid. Particles at these positions are displayed as 
small plus signs on the cathode ray tube and are then photographed. Thereafter, at each 
time step in the calculation, the particles are moved with velocities appropriate to their 
location and time, displayed on the cathode ray tube, and photographed. 

The distance a particle is moved is simply the product of the velocity at its current 
location and the time interval over which that velocity is assumed to exist (i. e. , the time 
step in the solution of the Navier-Stokes equations). Since the particles will not in gen- 
eral be located precisely at positions where the velocities are known from the solution of 
the Navier-Stokes equations, a way of estimating their velocities is needed. The com- 
ponents of the velocity of a particle are calculated as the weighted averages of the veloci- 
ties at the four closest positions at which those velocities are defined. Since u and v 
velocities are defined at different positions, velocities at eight positions will be involved 
in the movement of one particle. The weight assigned to a velocity is inversely propor- 
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tional to its distance from the particle in question. Figure 3 shows a typical particle 
and the velocities that are used to estimate the particle velocity. 

A film (C-271) entitled T ’Computer- Generated Flow-Visualization Motion Pictures” 
shows how this technique is used. A request card and a description of the film are in- 
cluded at the back of this report. 


CAVITY FLOW 
Background 

As a specific example of a problem that can be solved with the technique presented 
in this report, consider the startup. of flow in a long groove, over which an endless belt 
is drawn at constant velocity. A cutaway view of such a groove is shown in figure 4. 
This situation, which is difficult to study experimentally without disturbing the flow, has 
relevance in bearing and seal studies, where it is the limiting case of a spiral groove 
seal. In this case, the results of interest would be the integrated pressure over the 
belt, which is the net lift, and the minimum pressure, which would indicate whether 
cavitation would be a problem. 
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Figure 4. - Cutaway view of groove and belt. 


If the groove is long enough, the flow at each cross section, except near an end, 
will be the same. Steady flows such as this have been investigated experimentally and 
numerically so that comparison with prior work is possible. 

If the upper surface of a fluid-filled cavity is moving in its own plane with constant 
velocity, a circulatory motion of the fluid is set up within the cavity. Time-exposure 
photographs have been taken (refs. 8 and 9) of flows into which a tracer has been injected 
so that the qualitative features of the steady flow are known. For cavities with aspect 
ratios (i.e. , depth/width) of 1 or less, most of the fluid rotates about a point - the vor- 
tex center - where the vector velocity is zero. The main vortex occupies most of the 
cavity but small, weak, counterrotating vortices exist in both lower corners. For a 
cavity with an aspect ratio of 2, there is, in addition, a large, weak, counterrotating 
vortex occupying most of the lower portion of the cavity. 

Steady flow in a two-dimensional cavity has been analyzed (refs. 8 and 10 to 15) by 
numerically solving the steady Navier-Stokes equations. Burggraf (ref. 12) has also 
obtained an analytic solution to the linearized problem for an eddy bounded by a circular 
streamline. The starting flow problem, the fluid being initially at rest, was considered 
by Greenspan (ref. 16); however, only steady results were presented. 

The technique devised by Harlow (ref. 4) for calculating free surface flows was used 
in reference 5 to solve the problem of unsteady flow in a square cavity at a Reynolds 
number of 100. The calculations were carried out until velocities were no longer 
changing, at which point they were in excellent agreement with a numerical solution of 
the steady Navier-Stokes equations. In addition, the terminal position of the unsteady 
vortex center agreed well with the position of a vortex center estimated from a time- 
exposure photograph of a steady vortex. 
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The parameters of interest are the aspect ratio of the cavity and the Reynolds num- 
ber of the flow. Cavities with aspect ratios of 1/2, 1, and 2 are studied at a Reynolds 
number of 100. Flows with Reynolds numbers between 100 and 500 are investigated in a 
square cavity. 


Remarks 

For cavity flow it is convenient to choose the length and velocity of the moving wall 
as the reference length and velocity used to make the Navier-Stokes equations dimen- 
sionless. The cavity is assumed to be bounded by no-slip, impermeable walls. The 
fluid is at rest at the start of the calculation and the value for the initial uniform pres- 
sure is chosen to be unity. The reference location for pressure, where a value of unity 
is maintained, is the center of the wall opposite the moving wall. 

Space increments of 0.05 and a time step of 0.02 were found to be satisfactory for 
Reynolds numbers from 100 to 500. Some of the calculations were also performed with 
space and time increments halved, and essentially the same results were obtained. Ap- 
proximate stability criteria are discussed in appendix B. 

A difficulty with particle movement arose when a particle was near a wall where the 
boundary layer was very thin. It was noticed that particles tended to congregate along 
the upper part of the right wall. Then a crescent-shaped region devoid of particles 
formed along most of the remaining part of the right wall and extended into the cavity. 

In order to circumvent this problem the calculation of tangential velocity components of 
particles within half a mesh spacing of a wall was modified. 

A possible velocity profile in a thin boundary layer is marked "actual” in figure 5. 

If a particle is located at position P, for example, the interpolation scheme for particle 
movement using four mesh-point velocities underestimates particle velocity. If only the 
two mesh-point velocities within the cavity are used, particle velocity is overestimated. 
For the profiles shown in figure 5, however, the two-point interpolation is superior. 

This two-point interpolation scheme was used for particles within half a mesh spacing of 
a wall, and no further anomalous particle motion was noted. 

About 1/2 minute of IBM System/360 Model 67 computer time was required per di- 
mensionless time regardless of Reynolds number for a square cavity. However, longer 
runs were necessary to reach steady conditions at larger Reynolds numbers. The cri- 
terion of steady conditions we have used is that the position of the vortex center change 
less than 1 percent over a period of 5 dimensionless time units. For a Reynolds number 
of 100 this occurred at about a dimensionless time of 10; for a Reynolds number of 500, 
it was not reached until about a dimensionless time of 25. For cavities of aspect ratio 
1 '2 or 2. the number of mesh points in the longer direction was doubled and about twice 
as long a running time as for the square cavity was needed. 
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Figure 5. - Calculation of particle velocity in a thin boundary layer on a wall. 


Discussion of Results 

Cavity flow is characterized by the aspect ratio and the Reynolds number. Aspect 
ratio is the cavity depth divided by the cavity width. Reynolds number is the product of 
the length and velocity of the moving wall divided by the fluid viscosity. 

Aspect ratio, 1/2; Reynold s number^ 100 . - Figure 6 is a time-exposure photograph 
(ref. 8) of a steady vortex in a cavity with aspect ratio of 1/2 at a Reynolds number of 
100. Washing powder has been added to oil to make the flow visible. A tape is pulled 
across the cavity to provide the moving surface. Figure 7 compares the position of the 
vortex center estimated from this photograph with the position calculated by Zuk and 
Renkel (ref. 15) from the steady Navier-Stokes equations and with the position deter- 
mined from the unsteady Navier-Stokes equations at large time. The numerical solutions 
are in good agreement with the experimental value. The path of the instantaneous posi- 
tion of the vortex center is also shown; it starts under the midpoint of the moving wall, 
shifts downstream (i.e. , in the direction of movement of the wall) and into the cavity, 
and turns upstream slightly to attain its terminal position. 

Figures 8 and 9 compare terminal velocities from the solution of the unsteady 
Navier-Stokes equations with velocities calculated from the steady Navier-Stokes equa- 
tions (ref. 15). The comparisons are shown as velocity traverses through the vortex 
center; figure 8 shows velocity parallel to the moving wall in the vertical traverse and 
figure 9 shows velocity perpendicular to the moving wall in the horizontal traverse. In 
both figures the agreement is good. 
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VERTICAL POSITION 



Figure 6. - Time-exposure photograph of vortex in cavity of aspect ratio 1/2. Reynolds number, 100. (From ref. 8.) 


o EXPERIMENTALLY DETERMINED FROM FLOW 
VISUALIZATION STUDY (REF. 8) 

° NUMERICAL SOLUTION OF STEADY EQUA- 
TIONS (REF. 15) 

A NUMERICAL SOLUTION OF UNSTEADY 
EQUATIONS 

PATH OF INSTANTANEOUS POSITION OF 

VORTEX CENTER 


MOVING WALL 



HORIZONTAL POSITION 


O NUMERICAL SOLUTION OF STEADY EQUATIONS 
(REF. 15) 

NUMERICAL SOLUTION OF UNSTEADY EQUATIONS 

MOVING 



-0. 50 0 . 50 1. 00 

VELOCITY COMPONENT PARALLEL TO MOVING WALL 


Figure 7. - Position of vortex center in cavity of Figure 8. - Vertical velocity traverse through vortex center in 

aspect ratio 1/2. Reynolds number, 100. cavity of aspect ratio 1/2. Reynolds number, 100. 
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Figure 9. - Horizontal velocity traverse through vortex center in 
cavity of aspect ratio 1 12. Reynolds number, 100. 


MOVING WALL 



Figure 10. - Position of vortex 
Reynolds number, 100. 


START OF 
CALCULATION 


o NUMERICAL SOLUTION 
OF STEADY EQUATIONS 
(AUTHORS OF REF. 15) 
a NUMERICAL SOLUTION 
OF UNSTEADY EQUA- 
TIONS 

PATH OF INSTANTANEOUS 

POSITION OF VORTEX 
CENTER 


in cavity of aspect ratio 2. 


Aspect ratio, 2; Reynold s number, 100. - Two large vortices exist when the aspect 
ratio of the cavity is 2. Figure 10 compares the terminal positions of vortex centers 
calculated from solutions of the steady and unsteady Navier-Stokes equations at large 
times for a Reynolds number of 100. The authors of reference 15 calculated the steady 
case for comparison with the unsteady result. The numerical solutions are in good 
agreement with each other; no experimental value is available for comparison. The 
path of the instantaneous position of the upper vortex center is also shown and is similar 
to the path for cavities with aspect ratios of 1/2. 


16 




O NUMERICAL SOLUTION OF STEADY EQUATIONS 
(AUTHORS OF REF. 15) 


MOVING 

Vi/ALL 


LT) 

O 

Q_ 


QH 

LU 

> 


LuWEK 
vVmLL | 

50 


NUMERICAL SOLUTION OF UNSTEADY EQUATIONS 



VELOCITY COMPONENT PARALLEL TO MOVING WALL 


Figure 11. - Vertical velocity traverse through upper vortex 
center in cavity of aspect ratio 2. Reynolos number, 100. 



Figure 12. - Horizontal velocity traverse through upper vortex 
center in cavity of aspect ratio 2. Reynolds number, 100. 


Figures 11 to 14 compare traverses of terminal velocities calculated from the un- 
steady Navier-Stokes equations with velocities calculated from the steady Navier-Stokes 
equations. The agreement between the solutions is good. Figures 11 and 12 are tra- 
verses through the upper vortex center. It can be seen that the flow in the upper part of 
the cavity is similar to the flow in the cavity with an aspect ratio of 1/2 shown in figures 
8 and 9. Figures 13 and 14 are traverses through the lower vortex center. The flow in 
the lower part of the cavity is much slower (about two orders of magnitude) than in the 
upper part. 
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NUMERICAL SOLUTION OF STEADY EQUA- 
TIONS (AUTHORS OF REF. 15) 
NUMERICAL SOLUTION OF UNSTEADY 
EQUATIONS 



. 5 U • 

VELOCITY COMPONENT PARALLEL TO MOVING wALL xlO 2 


Figure 13. - Vertical velocity traverse through lower vortex 
center in cavity of aspect ratio 2 . Reynolds number, 100. 



Figure 14. - Horizontal velocity traverse through lower vortex 
center in cavity of aspect ratio 2. Reynolds number, 100. 



Figure 15. - Time-exposure photograph of vortex in square cavity. Reynolds number, 100. (From ref. 8.) 


Square cavity . - Figure 15 is a time-exposure photograph (ref. 8) of a steady vortex 
in a square cavity at a Reynolds number of 100 from which the position of the vortex cen- 
ter can be estimated. Figure 16 compares this estimate with the positions determined 
from two numerical solutions of the steady Navier-Stokes equations (refs. 12 and 15) and 
the unsteady Navier-Stokes equations at large times. The agreement among these vari- 
ous methods is excellent. The path of the instantaneous position of the vortex center is 
also shown; it is similar to the paths for cavities with aspect ratios of 1/2 and 2. 

Figures 17 and 18 compare terminal velocities from the solution of the unsteady 
Navier-Stokes equations with velocities calculated from the steady Navier-Stokes equa- 
tions (ref. 15). The curves, showing traverses through the vortex center, are similar 
to those already presented for a cavity with an aspect ratio of 1/2. The agreement be- 
tween the solutions is excellent. The flow in this square cavity is similar to the flows 
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VERTICAL POSITION 


o EXPERIMENTALLY DETERMINED FROM 

FLOW VISUALIZATION STUDY (REF. 8) 
o NUMERICAL SOLUTION OF STEADY 

EQUATIONS (REF. 15) 

□ NUMERICAL SOLUTION OF STEADY 

EQUATIONS (REF. 12) 

a NUMERICAL SOLUTION OF UNSTEADY 

EQUATIONS 

PATH OF INSTANTANEOUS POSITION 

OF VORTEX CENTER 


MOVING WALL 



Figure 16. - Position of vortex center in square cavity. Reynolds 
number, 100. 
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O NUMERICAL SOLUTION OF STEADY EQUATIONS 
(AUTHORS OF REF. 15/ 

— NUMERICAL SOLUTION OF UNSTEADY E^UAftONS 



VELOCITY COMPONENT PkRALLEL TO MOVING WALL 


Figure 17. - Vertical velocity traverse through vertex center 
in square cavity. Reynolds number, 100. 



o NUMERICAL SOLUTION OF STEADY EQUATIONS (REF. 15) 



Figure 18. - Horizontal velocity traverse through vortex center 
in square cavity. Reynolds number, 100. 
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Figure 19. - Vertical velocity traverse through vortex center 
in square cavity. 


in the cavity with an aspect ratio of 1/2 and the upper part of the cavity with an aspect 
ratio of 2. 

Figures 19 and 20 show how velocity traverses through the vortex center change as 
Reynolds number is increased from 200 to 400 in a square cavity. For larger Reynolds 
numbers the extent of the inviscid portion of the flow is greater. This is shown by the 
increased portion of the velocity profile that is linear about the vortex center. The only 
prior calculation available for comparison (ref. 12), a vertical traverse from a solution 
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Figure 20. - Horizontal velocity traverse through vortex center 
in square cavity. 


□ NUMERICAL SOLUTION OF STEADY EQUATIONS (REF. 12) 
NUMERICAL SOLUTION OF UNSTEADY EQUATIONS 



Figure 21. - Vertical velocity traverse through vortex center in 
square cavity. Reynolds number, 400. 


of the steady Navier-Stokes equations at a Reynolds number of 400, shows good agree- 
ment in figure 21 with the solution of the unsteady Navier-Stokes equations at large 
times. 

Figures 22 and 23 show velocity traverses through the vortex center in a square 
cavity for a Reynolds number of 500 at various times. It can be seen that at successive- 
ly larger times the inviscid portion of the flow occupies progressively larger fractions of 
the cavity. 
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Figure 22. - Vertical velocity traverses through vortex 
centers in square cavity. Reynolds number, 500. 



Figure 23. - Horizontal velocity traverses through vortex 
centers in square cavity. Reynolds number, 500. 














Figure 24 shows three-dimensional plots of the pressure fields for a Reynolds num- 
ber of 500 and the same times as the velocity traverses shown in figures 22 and 23. At 
the upper left corner of the cavity the pressure is below the reference pressure. The 
pressure rises along the moving wall, reaching a value greater than the reference pres- 
sure near the upper right corner on the last frame. A local minimum exists at the vor- 
tex center. The pressure fields for lower Reynolds numbers were similar to those 
shown in figure 24, but the mini mums and maximums along the moving wall were more 
pronounced. 

Figure 25 shows the effect of Reynolds number on the positions of the unsteady vor- 
tex centers as well as the terminal positions of the vortex centers for several Reynolds 
numbers from 100 to 500 in a square cavity. Increasing the Reynolds number shifts the 
terminal position of the vortex center upstream and farther into the cavity. These posi- 
tions fall quite close to the asymptote from linear theory (ref. 12), shown as the line 
drawn from the center of the cavity to the moving wall. The path of the instantaneous 
vortex centers moves farther downstream before turning down at the wall and then up- 
stream as Reynolds number increases. 

Flo w vi sualization. - Numerical flow visualization motion pictures were made for 
each of the runs. Selected frames from two motion pictures have been chosen to illus- 
trate the numerical flow visualization experiments. Figure 26 is for a square cavity at 
a Reynolds number of 100. The first frame shows the initial distribution of four parti- 
cles per cell at the start of the calculation. The remaining frames show that the fluid 
close to the upper moving surface is dragged along by the upper surface and forced to 
turn down at the right wall. The fluid flows along the right wall, gradually turning back 
to the left and then upward, only to be directly influenced by the moving surface and 
dragged to the right again. 

Frames from a positive print of the motion picture of flow in a cavity with an aspect 
ratio of 2 at a Reynolds number of 100 are given in figure 27. The first frame shows the 
initial distribution of one particle per cell. Succeeding frames show the development of 
the flow at dimensionless times of 1, 2, 3, 4, and 10. 

The circulatory nature of the flow is evident in these figures, even though it is seen 
much more dramatically in the motion pictures. The effect of the moving wall can be 
seen to propagate into the cavity with time. The lower portion of the cavity remains 
relatively quiescent, even at a dimensionless time of 10. Later frames show the begin- 
nings of the development of the weak vortices. The small triangle above the cavity is 
used to show the velocity of the moving wall in the motion picture. 
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Figure 25. - Path of instantaneous position of vortex centers in 
square cavity. 




(a) Time, 0. 



(c) Time, 2. 


Figure 26. - Numerical flow visualization experiment. 



(d) Time, 3. 

cavity; Reynolds number, 100, Dimensionless time. 
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(g) Time, 6. (h) Time, 7. 


Figure 26. - Continued. 
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(k) Time, 10. 
Figure 26. - Concluded. 
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(d) Time, 3. 


(e)Time, 4. 


(f) Time, 10, 


Figure 27. - Numerical flow visualization experiment Aspect ratio of cavity, 2; Reynolds number, 100. Dimensionless time. 






CONCLUDING REMARKS 


A computer program to solve the unsteady, two-dimensional, incompressible 
Navier-Stokes equations was written. The program was used to investigate unsteady 
flow in a rectangular cavity with a moving wall. The calculations started with the fluid 
at rest in the cavity and continued until no further changes in velocity occurred. 

Flow in cavities with aspect ratios of 1/2 and 1 and in the upper part of a cavity with 
an aspect ratio of 2 were similar at a Reynolds number of 100. This was true of the un- 
steady, as well as the steady, flow. 

As the Reynolds number increased from 100 to 500 in a square cavity the inviscid 
portion of the flow occupied a progressively larger portion of the cavity. At all these 
Reynolds numbers the inviscid portion of the flow about the vortex center increased 
during the initial period. 

For a square cavity the path of the instantaneous position of the vortex center ap- 
proached the downstream wall more closely with increasing Reynolds number, even 
though the terminal position of the vortex center moved closer to the upstream wall. 

No stability problems were encountered in any of the runs. The velocities at large 
times were compared where possible to velocities calculated from the steady Navier- 
Stokes equations or the results of steady experiment; generally satisfactory agreement 
was shown. 

A numerical flow visualization experiment was described. This technique uses 
marked particles to show the motion of the fluid. A microfilm recorder is required in 
order to display and photograph the particles. At the end of a run the sequence of photo- 
graphs can be projected as a motion picture which shows clearly the development of the 
flow. The coding necessary to perform such an experiment is included in the computer 
program for solving the Navier-Stokes equations. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, February 1, 1971, 

129-01. 
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APPENDIX A 


SYMBOLS 

Mathematical FORTRAN Description 

symbol symbol 

d D correction term, eq. (13) 

g v , g body force in x- and y-directions 

x y 

g x , gy GX, GY dimensionless body forces, Lg/ 

h H dimensionless cavity length in y-direction 

i I subscript indicating y-position 

j J subscript indicating x-position 

L reference length (length of moving wall for cavity flow) 

A. 1 _ 

n N superscript indicating n in time step 

p pressure 

J.U 

q superscript indicating q n iteration 

r R abbreviation, eq. (17) 

t time 

t T dimensionless time, Wt/L 

u velocity in x-direction 

u U dimensionless velocity in x-direction, u f W 

v velocity in y-direction 

v V dimensionless velocity in y-direction, v/w 

W reference velocity (velocity of moving wall for cavity 

flow) 

w W dimensionless cavity length in x-direction 

x horizontal direction 

x X dimensionless horizontal direction, x/ L 

y vertical direction 

y Y dimensionless vertical direction, y j L 

6t DELT time step 
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Mathematical 

symbol 

FORTRAN 

symbol 

Description 

6x 

DELX 

space increment in x-direction 

6y 

DELY 

space increment in y-direction 

~V 


fluid kinematic viscosity 

V 

vise 

dimensionless kinematic viscosity 

p 


fluid density 

<p 

PHI 

dimensionless pressure, p/pW^ 


OMEGA 

relaxation factor 


t 


f 
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APPENDIX B 


DESCRIPTION OF COMPUTER PROGRAM 


The three major calculations performed by the computer program are the solutions 
of the equations for pressure and for both velocities. Additional calculations are re- 
quired to determine the position of the vortex center and for particle movement, plotting 
routines, etc. Since the continuity equation will not be satisfied exactly, the d- . will 
not in general be identically zero. It is necessary, however, that the continuity equation 
be satisfied sufficiently closely at each time step. Following Harlow (ref. 4), we make 

_3 

sure that each d. . is smaller in absolute value than 3.5x10 
J 

The Poisson's equation for pressure will be solved iteratively. In order to expedite 
convergence, successive overrelaxation (see ref. 17) is used. With the introduction of a 
relaxation factor co, equation (16) becomes 





(Bl) 


where the superscript q indicates the q iteration. The most recently computed 
values of the <p shown without superscripts are used. A relaxation factor of 1. 8 has 
been used without any attempt at optimization. 

The iteration is continued until the field converges to some desired degree. Satis- 
factory results have been obtained when the criterion suggested by Harlow (ref. 4) 


<Pi 


? . - 1 


M fill 


< 2x10 


-4 






<1 


2 2 
+ U i,j + V i,j + 


(B2) 


|g x w 


lg/1 


is satisfied at each point. 

The velocities are calculated explicitly so that the time step is limited by stability 
conditions. Hirt (ref. 18) has suggested approximate criteria 


and 


v > — 6t u 2 
2 


(B3) 


x 1 _ 2 3u 
v > — <5x — 

2 3x 


(B4) 
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where u is the average maximum fluid speed and Ou/?x is the average maximum ve- 
locity gradient in the direction of flow. Equation (B3) indicates a time step of 0.02 for 
a Reynolds number of 100. Equal space increments were used to give a 20x20 mesh and 
it was observed that equation (B4) was satisfied also. The same time and space incre- 
ments were also found to be satisfactory for Reynolds numbers as large as 500, indi- 
cating that these criteria are somewhat conservative for this cavity flow. 

It remains to expand terms such as u^ j prior to coding. Capital letters are used 
to represent the FORTRAN notation for the finite difference variables. Since fractional 
subscripts are not allowed, U(I, J) is used for Uj j + \/2 arM * V(I, J) is used for v j + j/ 2 j> 
the other variables are cell-centered, so that arguments (I, J) represent (i, j). 

The computer program was written to be run on an IBM 360 Model 67 computer in 
the time-sharing mode. Unless the program is being restarted, the BLOCK DATA sub- 
program named CBLKC must be loaded before the main program named MAIN. Addi- 
tional routines are required for calculating running time (CPUTIM), microfilm plotting 
(LRGON, LRCNVT, LRGRID, L RANGE, LRCHSZ, LRLEGN, LRXLEG, LRYLEG, 
LRTLEG, LRCURV, LRMRGN, LRMON, LRMOFF, LRNON, LRNOFF, LRSON, 
LRSOFF) given in reference 19, and drawing isometric views (PLOT3D, ROTATE) 
given in reference 20. 


Outline of Computing Procedure 


An outline of MAIN is given below. The names of SUBROUTINE or ENTRY calls 
are listed at the right. 


(1) Read and write the input data. 

(2) Calculate constants that will be needed. 

(3) Calculate and plot the initial particle positions. 

(4) Calculate D and R at time T. 

(5) Calculate PHI at time T + DELT. 

(6) Calculate U and V at time T + DELT. 

(7) Calculate the position of the vortex center. 

(8) Calculate and plot the new particle positions. 

(9) Print D, U, V, and PHI if T = 0, 1, 2, 3, . . . 

(10) If T < TLAST, return to step 4; if T = TLAST, 

prepare final prints and plots. 


INIT, CONSTT 
PART, PRINTI 
DANDR 
PRESS, ITER 
UANDV 
VTXCTR 
MOVE 
PRNTN 
PRNTF 
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Input 


NAMELIST is used for data input; the variables, with default values listed at the 
right, are any or all of the following: 

DELX, DELY 

space increments 

0.05 

DELT 

time step 

0.02 

W 

depth/width ratio of cavity 

1.0 

TLAST 

dimensionless time at which to end calculations 

2.0 

vise 

kinematic viscosity 

0.01 

OMEGA 

relaxation factor 

1.8 

MOVIE 

logical variable, true if a motion picture is desired 

.FALSE. 

SECOND 

logical variable, true if program is being restarted 

.FALSE. 

NPPC 

number of marked particles per cell 

1 


Output 

Marked-particle positions are recorded at every time step in the calculation. 
Printed output, consisting of all values of D, U, V, and PHI is generated at T = 0, 1, 

2, 3, . . . . At the end of the calculations, microfilm plots are made of the instanta- 
neous positions of the vortex center. Horizontal and vertical velocity traverses through 
the vortex center and an isometric view of the pressure field are also plotted. 

The printed output from the final time for the default values of the input variables is 
given in appendix C. 
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APPENDIX C 


OUTPUT FROM SAMPLE PROBLEM AT DIMENSIONLESS TIME 2 


-p. 21A6E-93 

-0. 155AF-C3 

-0. 7725F-CA 

-C.l 55AF-33 

-p * A 3 R 7 F — 9 A 

. A P 5 P E- ' A 

r . 4 3 ? 91 -54 

: . P46i-'. 3 

". 2A7RT-U 3 

362PF-E 3 

9. ip f. AC— 0 3 

C.2 7R9F-03 

P..9A6 3R-04 

C . 7 62 IF- 04 

C . 77A? c — OA 

C. 940 3 r- JA 

-P. 11 4 4i= -:a 

-p.4 7 7 I f — . A 

-0 .74 39 r-JA 

-r) . 4 >»7 1F-0A 

-9. 7 7 ? 5f-C A 

C.163RE-OA 

-0.6294F-04 

'3 .5531 F — 94 

.6967F-CA 

C.l (‘RIF- C 3 

J . >’P/,9E- JA 

".11 7 C F - C 3 

J • l 197F-U ^ 

J.22AAF-L 3 

f- . ^7 A Af -0 A 

-C.7A7AE-0A 

-C. 454PE-04 

0. t 319F-03 

C . 1 90 6 £ 3 

9. * 141=- 34 

-(j . r 4 r r r: - : a 

p . ll aat - i <* 

. .915SF-',4 

L. 47141'-.' A 

- J . 9 3*6 F— 0 A 

-0.3 77* 7 E -04 

-C .l?16F-t)A 

C . 13A3F-0A 

r, 1161F-03 

C. Ap«-«7r-_A 

", 0 7 751'- "A 

>. 12 76F-C 3 

' .?2 15F-C t 

■".4 35 7 r -0 4 

o. 3 1 6 5F — 0 A 

0.2956E-94 

C.1155E-03 

P. 1 IA9C-G3 

C . i 1 7* 2 F - 3 3 

-C . 1 5 7AE-03 

-:. 1CA9F-3A 

- ' .2671F-C4 

... 221 /E-j a 

P. 3 6 4 R f — C 4 

-". 713AE-04 

-C.511AE-04 

C. 152CF-04 

C.079PC-OA 

-C . 3165C-04 

C . 7 " 9 9 f - 34 

g . i 

■9.11 7 = 1 -C ? 

C . ?C 2 l F- 3 

P.91 73r-C A 

-r. R73F-0A 

0.2009E-OA 

0.5347F-04 

0.5^11 E-DA 

-P .39 jAF-C A 

- C • ?2 P 3=- 53 

-(_ . 69 6 2r — 3 A 

3. 2 3 = 4 E-P A 

9.14 31 t-J> 

-J,7A°0F— C A 

-C. 1A22F-03 

0.5966E-94 

9. 43 1 6P-0 A 

0.1 31 IE-C4 

r .1 IR0E-C3 

; .6 L74F- 34 

G . 3C99F -;A 

0 . 7 79 4f-L 3 

3. 72 4.H- jA 

9. 3 2 7HF -C 5 

-p .421 AF-9A 

0 .7A2RF-0A 

C.2? 6 6F -05 

0. A 2 2 6F - 0 A 

— C . 7 59 3F -0 3 

-C . L 779C-03 

P.536AF-1A 

3 . 1 7 p p F -c A 

. ? 5 7 5 ”— J 4 

-■.95 37r-,,5 

-3. 3320F-CA 

C.34H7E-CA 

-C .7 1C^F-0A 

0 . 7RP?r-GA 

C.l 1 9 2F-P 5 

0 • a P 2 it — CA 

C. 1633F- >3 

j • 6APAF-CA 

D.R224F-0-* 

0.1 550 1 —F A 

0. A?«2E-0A 

C.2CC9F-CA 

-0. 6?88r-0A 

-0.1 160F-93 

-C • 2 B 2 5 C -« 3 

-C .9437E-3 5 

" • 33 36 r — jA 

'.?9£-cr-r4- 

'.11 4 A F — 0 A 

9. 67. 74F-0R 

3 . 1651F-0A 

-0.1 101F-03 

0. 8154=-C4 

-C. 1 40 l*=- 04 

C. 1 1 P9F-L 3 

C. 61 26=-J4 

C. 61 76F-.A 

.-.446AC-' 4 

3 • 4 2 44 f— U A 

-0. 1 99 7 f -0 A 

-O.23PAE-06 

-9.67 4 IE -0 4 

-C.935RR-9A 

-0.1 733E-J3 

-C • 2 3R AT-C 3 

C .U6PE-0 3 

9. 1 7MRP-3A 

. A 1 9 7 r -04 

-' . Ida o E-i'4 

-0 • ? 4 a 2 r — f 5 

-0. 3827E-0A 

0.414 3E- 94 

-0. 1532F-0A 

9.4 3 94 c - PA 

-C.6795E-9A 

C .8 34 5F-0A 

C* .92 3°F — 14 

9 . A6A9r- r u 

?.? 3 3 A F - 1, A 

C. 1 l 9' , F - C A 

-O.A107F-CA 

-0.120PE-03 

-0.11P0F-03 

-P.ARPPl-OA 

-C-. 1 66°F -0 3 

P. IA31F-D3 

-0. ? 26 5F - " A 

-9. 1 527 E -1 A 

9. j 

-P ,44«3F-f 6 

-0. 213AE-0A 

-0. 17RPF-04 

C.3RHPF-0A 

-0 .3451 F-OA 

0 . ? 319E-jA 

C . 7 39 1 F- CA 

P. 1 907F - j A 

p. 4 77 FT -04 

J. 5 040 F- j 6 

-P . 46 76 r — L A 

-0 . 2 68 R F-0 A 

-0 • 1 23 3E-03 

-0.53A7 r -0A 

-9. 9 D6CC-0A 

-0. 3 4 5 7 F —0 A 

U.1776E-C3 

f .AC 4 3 F ->4 

-3 . 1 33«-r-C4 

~ . 1007F-V.5 

-r. 1 A 3 1 F- C 6 

-0. 362AF-0A 

0.3660E-04 

-0 .2217P-0A 

- n . 5 3 3 5 l -04 

C . PA9AF-0A 

-( .3PL5r-C4 

p. a?97f-:a 

0.0 4 3 7Z~k. 5 

G • “ 3 A 5 f-C4 

c . 15 37F-5A 


-9. 1 6*3 3F-04 
0. 1 142F-04 
-0. 7471F-04 
-r . 3PA4F-94 
-C. 7737E-04 
-0. '»1A?F-0A 
-n . lr 33E-03 
-i . f o a a r -o a 

-0. 1 A ^ 9F — OA 
-0. 6I62F-34 
-0. 1 3 70E-0A 
-0. 46 4 AT -04 
-5 .7.126F-04 

- J. 3fi7ir-04 
-0. Af OC F - 04 

- >, 5797E-"4 
-0. HALF -04 

3« 75F-04 
7 .'<49 3 F-0 5 
-3. 42* 4= -04 
-r. ?3l2f-94 


•C. 1499E-03 
-P.510BE-94 
•0.6717E-34 
-0. 7420F-C4 
-0.91G7E-94 
■0. 307&E-04 
C.2980F-06 
■0 .2L°7F-D4 
-C.93RRF-94 
0.17? 3E-04 
-0.7AH3E-0A 
C. 3754E-04 
■0.3 31 (SE-OA 
■0 . I 9ft 7 F -0 A 
-G.4L47E-04 
C .5960 E-3 6 
*■''.4711 E— 0 A 
■ C . 74R5F -04 
■0.6 40 7 F - 3 A 
-0. 1 3526-04 
-r. 3 1 3 56 - 04 


■C.S75CF-0A -n.S245F-04 
■0.5 46 Ch -0 A 0. 16756-04 


■0 • 466 A t;— 0 A 
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APPENDIX D 


COMPUTER PROGRAM LISTING 


SflUTt'V. ' F i)\STF»,)V NAVl ru -S l^fFS ’ 1. 1 , A I I " ► < 
P' " AV I r Y M ".‘1 


r cnn nr .r.qjfr.r 

c o;Fr f v i v , y p , o s a- A w f = o a v it/ , > i o r | ■ ' n = l H 1 L ! ‘* 

c t_ r* a i > r ’ i kc , r *m t N 

r ijsf '*w . k ruiMir.rn')[j*n 

r p»im . pin k» : ct , , , r hit , F‘>ast 

r m a sf >*. sr \ r > t. r no i h ci 

r 

r sf.i-! If p 'K f Al r» »| a 1 | n*>! 


C 

r 

r-- — inf w[CHf’Fir a ‘ j t r plot pthfcts nil hcav*'** i -* f s • r n<i <-axis 

r j r mi ^ if k r, F i u « f' tr put is SHf «■*. as a m :*■ l * . a • ’ *■ .• • • < r -“ r l ■*. " 

l nr.fr A l “.TV! F , sr C TM’ 

C*" n f i x , nf;i v , ori_ t, ^ , h, ti a c t, v! sr , * o \ ,r. x . .v , 

1 . [ L*S 1 ,Of L v SO. • WA ,GKOHK , CVOf L V, VI Sr < , v l^ r V ,f i *"»r , vr..i ?, 

? Xi > v, n vs, V X , V Y , V* Y ,Hif I X ,» Of t Y , T = ' *0 Y, ) , ■ )V IX , Sr v ,M AST , 

T V ,L , II f J.J , I I PI , JJ°! , 1 IM I t . 1 JM I ,M # I, >£ I . 1 - 1 1 ,-»w\v ,r |-f «-V X, X" AX , 'IP A>.* * 

C-' w"riJ/ljVP/U ( , AA ,? ) , Vf A*. » AA , H , PM ! ( AA , A A | ( '♦'< ,'|A I • * ( A '* . A A 1 , 

1 |>f Of , PC PN, '-‘AX I T, 1 Tf RP, I ir« | , S r*Ax 

Cr MW l'K / m ]/•*, ML AST , X ( lis u r ) 

C'"* ,J ic m/vc/v/p x ( i *500 » , vr v * i ‘i mi , j-mc s » i , 

Cl' ^ / PL /VPHTFL I A ) ) , i |p hr* I c sf I , VF •< r PS < ■> 1 I ,l* °F n< 'l t ii I 

\A yri IST/INPlJT/riU X , rir l Y ,'iri T , «, Tl AST , V I V.,f "C ‘SA , l 'f V It , SH M\0,',P>»'‘ 

TAIL f.P'ITl *1 Ml) 


f Vise IS T lit Kter 1 ' HM PF *MF a F^ MRS Ml’ 14 ! 

r 1 ' v F S A is TH[. flfi. ax M Ii'\ FArrnP 

f N-’PC IS nr \IJVAI-* ( F P A » F T C l F S J F ■* CHI 

r srrn K( ) is .f»uf. ir phiopam js Aifins m sta* i 


!^ THIS CAST, n f ■ • T IPftF «Vl f'C x TATA 
L <- r . CS 
OF l V =0 .C ^ 

Oi- l T = 0.0? 

« = 1 . C 
H* J .( 

T| A S T = ? . 0 

vr so=o. oi 

i j v rS A = t ,P 
VI'VI r =.F ALSt . 

Sf n= . F Al. SF. 

N J P r " 1 
As 1. C 
AC *"■'"» = A 

— k Ail A NO Wn! ITT TH 4 ' OAT A 
| N p I J T 1 

XK IT* (0,10 ) OF I X, -jrLV.OEl T f w,H,T LAST, vise, r-^rSA, OX, r,Y 

ifisfcondi ^tax I'.i r ,i», v',p«t ,n , l , k , nr.f'N, Pd'N, max i r , n fpo, r r, - 

1 MlAST,MMAx,y,Y t AC Of, , VF P T p S , VO 7 'F l ,H -'«< J (»S , 1 1 P S’ Orj , V ( X , V Y 
IF (^F CONTI OF MI NO K 


I V-T , YsO) I X- W 


l F p r sAII ■< V s 


( X O , Y = H I ( X = W,V=M | 



o n 


Ifi .MOT. SECOND) T=0.0 

C CALCULATE THE CONSTANTS NEEDED IN THE VARIOUS SUBROUTINES 

CALL CONSTT 
IF i SECOND) GO TO 25 
TF IRST = 0 .0 

PRINT IDENTIFYING INFORMATION ABOUT THE RUN AND DESCRIPTIVE 

MATERIAL FOR THF MOVIE 

5 CALL l N I T 
IF(MOVIF) CALL LEADER 
I F { SECOND) GO TO 6 

C CALCULATE INITIAL PARTICLE POSITIONS 

CALL PART 

C PRINT INITIAL PARTICLE POSITIONS 

CAtL PRNTI 

C START OF CALCULATICN LOOP 

L = I 
K = 2 

TF I R S T=DEL T 
NF IR ST= 1 
GO TO 7 

6 TF I R ST= T + DE L T 

NF I R ST= ( T+DELT+O. 00 I ) /DELT 

7 DO 14 N^NFIRSTtNLAST 

T- TF IRST+FLOAT ( N-NF IRST ) #DFLT 

C CALCULATE D AND R FROM U AND V AT TIME T 

CALL OANDR 

C CALCULATE PHI AT TIME T+DELT 

CALL PRFSS 

C CALCULATE U AND V AT TIME T +OFLT 

CALL UANDV 

C CALCULATE THE POSITION OF THF VORTEX CENTER 

CALL VTXCTR 

C MOVE PARTICLFS 

CALL MOV F 

C PRINT NEW PARTICLE POSITIONS 

CALL PRNTN 

8 IF ( ABS { T-ACHG) .LT. 0.001) GO TO 17 

13 GO TO ( 1 I * 1 2 ) * L 

11 L = 2 
K* 1 

GO TO 14 

12 L = 1 
K - 2 

14 CO NT I NU F 
GO TO 9 

17 CALL CPUTIMIN2) 

C Rtj NT I M IS CPU TIME IN MINUTES 

RUNTIM=FL0AT(N2-N1) /lOOOCOO.O 
WR I T E ( 6 t 15 ) T 
WRITE(6,28) RUNTIM 

WRITE ( 6 v 2 3 ) ((D(I,J) f I=3,JJ),J=3,JJ) 

WR ITF<6 , 1 8 ) UU ( I * -J«K ) « J=*,JJ ) v I = 3« II) 

WRITE <6, 19) (( V( I t .I,K) ? J=3 t JJ ) f I =3 t I I ) 

WRITE (6, ?0) ( ( PH I ( I f J),J=?,JJPl) ,1=2, I I PI) 

WR IT E ( 6 t 26 ) 

BACKSPACE 6 

achg=achg+a‘ 

WRIThdul r,U»V»PHI,O f LtKf DCO N , P C ON , M A X I T , I TERO , I TFR I , - 
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I MLAST.MMAX, X, Y.ACHG, VERTPS , VPR0FL,H0R7PS ,UPROFL ,VCX,VCY 
REWIND LC 

I F ( T.LT.TLAST ) GO TC 13 
9 IF(MCIVIE) call tail 

WR ITE( 10 ) T,U, V,PHI,D,L,K,DCON,PCON,M4XI T, ( TFR 0*1 TER [ , 

1 MLASTt MPAX,X,Y, ACHG, VE RT PS , VPROFl , HORZPS , UP RO FL ,VCX,VCY 
REWIND 10 

C PRINT FINAL FRAME OF MOVIF ANO POSITIONS OF VORTFX CENTERS 

CALL PRNTF 
STOP 

25 WR I TE(6, 15 > T 

WR ITE(6*23) (IDl l , J ) , J=3 » J J ) , 1 = 3, II) * 

WRITE! 6, IB) ( ( U ( I » J » K ) , J = 3 , J J ) ,1=3,11) 

WR ITF(6, 19 ) ! I V ( I,J,K),J=3,JJ),I=3,II) 

WRITE (6, 20) ( (PHI ( I,J) , J=7, JJP1), 1=2, II PI ) 

WR ITF(6,26) 

BACKSPACE 6 
GO TO 5 

15 FORMAT! 4HIT = ,F10.4I 

16 FORMAT ( IHO, 10(F6. 3) ) 

IB F0RMAT(2H1U/ (10E12.4) ) 

19 FORMAT ( 2H1V/C 10E12.4) > 

20 FORMAT (4Hl°HI/ (llEll.4) ) 

23 FORMAT ( 2H0D / (10E12.4)) 

26 FORMAT ( 1 X, • JOB TERMINATED HERE*) 

28 FORMAT (12H0CPU TIME = , E 1 6 . B , 7HM INU TE S ) 

END 


FUNCTION FPR FSS(X,Y) 

COMMON DFLX, DELY, DELT, W, H, TL AST , V I SC , OMEG A , GX, GY , 

1 OELXSO, DELY SQ»DX4,DY^»GX DFLX , GY06LY » V I SCX »V ISCY ,COMOG, OMGO 7 , 

2 XOY, XOYS, VX,VY, VX Y, FDE L X, FOE L Y.TDXDY, 0X00 Y,DYODX,GYHGXW,Nl AST, 

3 K,L,II,JJ,IIP1,JJPI,I I M 1 , J JM 1 , N, T , DEL ph I , DM AX , CPNMA X, XM AX , NPAR T 
COMMON/UVP/U! 44 ,44 , 2 ) , V ( 44, 44 , 2 ) , PH I (44 ,44 ) , 0! 44 ,44 ) , P ( 44, 44 ) , 

1 DCON, PC ON » M AX I T » I TERO, I TER I , SO«AX 
I = Y 
J = X 

FPRCSS=5C.0*PHt ( I, J) 

RETURN 

END 


SUBROOT INF CONSTT 

C CALCULATES CONSTANTS 

COMMON OELX , DELY, DEL T , W ,H , TL AST , V I SC , OMEG A , GX , GY , 

1 DFLXSO,DELYSO,DX4,DY4,GXDELX,GYOELY,VISCX,VISCY,COMOG,OMG07, 

2 XOY , XOYS, VX » VY »VXY »FDEt X , FDELY , T OX DY, DXOOY , DYOOX, GYHGXW , NL A ST , 

3 K ,L, II , JJ, II PI, JJP1,I I Ml, JJM1 ,N,T, DELPHI »DMAX *CONMAX »XMAX , NPART 
COMMON/UVP/U (44,44,2), V! 44, 44, 2) »PH I ( 44 , 44 ) , 0( 44 ,44) ,R ( 44 , 44) , 

l DCON, PC ON, MAX IT, I TERO, I TER I , SOM AX 
GX=0 . 0 
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GY =0 .0 

DELXSQ=DELX**2 
DELYSC=DEL Y**2 
0X4=4. 0*DELXSQ 
DY4=4.0*DELYSQ 
GX DEL X= GX#DEL X 

gydely=gy*dely 

VI SC X=2 . 0* VI SC/DELX 
V I SC Y= 2 .0*V I SC/DEl. Y 

C OMEGA IS THF RELAXATION FACTOR 

COMOG- 1.0 - OMEGA 

0MG0Z=0MEGA/(2.0*( l .O/DEL XSQ+ 1 . 0 /DEL YSQ ) » 
X0Y=DELX/(2 ,C*DELY ) 

XOYS=VISC*DELX/DELYSQ 
VY=V I SC/ DELYSO 
VX=VI SC/DELXSQ 
VXY=VI SC / ( DELX*DEl. Y) 

FDELX=4.0*DELX 

FDELY=4.0*DELY 

TDXDY= 2 . 0*DEL X*DEL Y 

DXODY=DELX/DELY 

0Y0DX = DELY/DELX 

GYHGXW=ARSIGY*H )+ARS(GX*WJ 

NLAST=(TLAST +0.00 I I/DELT 

11= 2.0 t H/OELY 

JJ= 2.0 + W/DELX 

I I PI = I I +1 

I I Ml — f t — 1 

JJP1=JJ+1 

JJV1=JJ-I 

RETURN 

END 


SURROUTINE INIT 

C PRINTS INDENTIFYING INFORMATION 

INTEGER BBBR ( 8 ) 

EXTERNAL FPRESS 
LOGICAl CUBE 

COMMON DELX,DELY,DELT, W,H,TLAST, VTSC,OMEGA,GX,GY, 

1 DELXSOtDELYSQf DX4, DY4, GXDEL X, GYDEL Y, VI SC X, VI SC Y ,C OMOG ,OMGOZ, 

2 XOY.XOYS, VX ,VY,VXY,FDELX f FOEl Y.TDXDY, DXODY.DYOOX, GYHGXW.NL AST, 

3 K,L,II,JJ,IIPI,JJP1,IIM1,JJM1,N,T,DELPHI,DMAX ,CONMAX,XMAX , NPART 
COMMON/ UVP/U (44,44,2 ) ,V< 44,44, 2 J ,PHI( 44,44) ,D( 44,44) , R ( 44,441 , 

I DC ON, PC CN, MAXI T »ITERO» ITERI »SDMAX 
COMMON/M I /M, ML A ST, X( 1600),Y( 1600) 

COMMON/ VC/ VCX (1500 ) ,VCY( 1500 ) , JCH< 50 ),NPPC 
COMMON/LABEL /LABI 15) 

COMMON/ PL/VPR0FL(50) ,UPROFL ( 5C ) , VE RTPS ( 50 ) ,HORZPS( 50) 

DIMENSION LARAU5) 

DIMENSION S ( 3800 ) 

DIMENSION CHART ( 3 ) 

DIMENSION FI N(6 ) 

DIMENSION VXC( 16) 

DIMFNS ION BX ( 5 ) , BY ( 5 ) 

OTMFNSION DE(5) , CE (6 ) , I 0 ( 3 ) , IK (3 ) 
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I 

I 


| DIMENSION INU ( 7) ,TNT2C 20) » INT3( 18) , INT4( 14) ,UNX (?) .(JNY (?) 

I DIMENSION T ITLl ( 30) 

| DIMFNSION FRAMX(8) ,FRAMY(8 ) 

I DIMENSION XV EL ( 1 1 ) , YPOS< 5) ,TPAR( 12) * YVE L ( L2>, X POS( 5), TPN 0(13) 

I DIMENSION VCXP(t600),VCYP< 1600) 

,1 CURE =. FALSE. 

j DATA LARA/4H PR E » 4HS SUR , 4HE . , 4HWALL »4H MCV,4HES l,4HN Y ,4H0IRF, 

|. 1 4HC T 1 0 » 4 HN . / 

t DATA FRAMX/1.0,9.0,1.0,9.9,1.C,1.0,9.0,9.0/ 

| DATA FR AMY/ 8. 066, 8.066, 1.934, 1 .934, 1.934, 8.066,1.934,8.066/ 

| DATA VXC/’ POSITION OF VORTFX CENTERS 1UPPFR WALL MOVES FROM LEFT T 

I 10 RIGHT'/ 

! DATA FIN/' $C9$D4£R9$R9$R6THF F N D ' / 

I DATA SYMBOL/'*-'/ 

I DATA WL / ' 3 ' / 

f DATA CHARX/'X'/ 

I DATA CHARY/'Y'/ 

t: DATA CHART( 1 ) /• T= • / 

DATA RX (1) ,BX(4) ,BX( 5 J /3*2.5/,BY(l) ,BY{ 2) ,8Y( 5>/3*2. 5/ 

,i DATA XW/2.0/,YW/2.5/ 

!l DATA DF/'SQMAX = '/,CE/« CONMAX = •/ 

f' DATA IO/'ITERO = '/,IK/'ITERI = •/ 

DATA INTI/' £03tR9$R9IR4$S3I NPUT OATA'/ 

DATA INT2/' $D5 *R 7 *GND$GF X = tL OXXX XXXXXTR6$GNDtC,FY =tKl XXXXXXXXTR6 - 
! l $GND$GFT = SL 3 ' / 

| DATA INT3/' $D7$R7R?WNFtWF = 1/VISCOSITY = tL6 t R40M EGA = - 

1 $L 6 ' / 

DATA INT4/' 5D9 $R7G tWNX f W F = $L 5XXXX $R 6GTWN Y$WF = $L7'/ 

DATA Tl Tl.1/' *C4t.R9tR9$R9A$03t t 8C0MPUTER SOLUTI CN$D3$L9 tL20F THE $C3 
ltR9$P. 3UNSTFADY N AV IER- S TOK E S FQOAT [0NS$D3tL9$L9$L?F OR FLOW*D3$L6IN 

2 A$D3SL9$L4TW0 DIMENSIONAL CAVITY'/ 

DATA UNX/22. 0,32.0/, UN Y/48. 0,48.0/ 

DATA SVCU/'H'/ 

DATA XVE L / ' VELOCITY COMPONENT PARALLEL TO MOVING WAIL’/ 

DATA YPOS/ ' VERT ICAL POSITION'/ 

DATA TPAP/* VERTICAL VELOCITY TRAVFRSF THROUGH VORT Ex CENTER'/ 

DATA YVFL/ ' VELOC 1 TY COMPONENT PERPENDICULAR TO MOVING WALL'/ 

DATA XPOS/ • HOR I 7.0NTAL POSITION'/ 

DATA TPND/'HORI ZONTAL VELOCITY TRAVERS THRnuGH Vf’RT r < (ENTER'/ 

DATA BBRR/5* ' ',' LEO',' DON','OVAN'/ 

VPROFL (l)=l .0 
HORZ PS ( 11=0.0 
VERTPS ( 1 )=W 
DO l 1=1,15 
1 LAB ( I ) =L ABA ( I ) 

IRF= 1 .O/VI SC 
I I M2 = 11 -2 
J J M 2= J J-2 

- FIIM?=IIM2 

FJ JM2 = JJM2 
FI l =1 I 

„ FJ J= J J 

XRNG=JJ+3 
YRNG=l 1+3 
IF(W-H) 6,7,8 
6 XMRG=3.5 
YMRG=2. 0 
GO TO 9 
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7 X M R G =2 . 0 

ymrg=2. 0 

GO TO <7 
a X MR G=2 . 0 
YMRG=3. 5 

9 RX(2)=FLOAT(JJ )+0.5 
BXm = 8X(2 » 

RY( 3)=FLT1AT( Il)+0. 5 
BY ( 4 ) = BY ( 3 ) 

— PRINTS IDENTIFYING INFORMATION 


CALL l RGON 

CALL L RCNVT I DFL X » 3 » I NT 2 (6) ,3*8,2 ) 

CALL LPCNVTIDELY,3,INT2 (12) ,3,3,2) 

CALL LPCNVTI PELT, 3, INT2I 19) ,3,5,2) 

CALL LRCNVTI I RE , l , I NT 3 I 9 ) , 1 , ft , 0 I 
CALL LPCNVTI OMEGA* 3 , I NT3 ( I'M ,3,5.2) 

CALL LRCNVTI GX, 3, I NT 41 7), 3, 4, 2) 

C A I L LRCNVTI GY, 3, I NT 4(14), 3, A,?) 

CALL LRGRIDI I, 1,0. 0,0.0) 

CALL IRA NGEIC.O, 55.0,0. 0,55.0) 

CALL LR LEG N( INTI, 25, 0,1. 0,9. 6, 0.0) 

CALL LRC URV( UN X, UNY, 2, 2, SYM.O. 0) 

CALL LRLEGNI I NT 2, 33,0, I. 0,9. 6, 0.0) 

CALL LRL EG Nil NT), 72, 0,1. 0,9. 6, 0.0) 

CALL LRLFGN UNT4,55,0,l.C,9.fc,L.0) 

CALL LRLEGNI T ITL 1,153,0,1.0,9.6,1.3) 

RETURN 

ENTRY LEADER 
CALL LRMON 

CALL LRMRGNI 0. 0 ,C. 0, O.C ,0. 0) 

CALL LRANGEIO.C, 1C. 0,0. 0,10.0) 

CALL LRGRIDI 1 , 1 ,0 .9 ,0.0 ) 

CALL L R NON 
DO ?C 1=1,32 

20 CALL LRCURVI X,Y ,0, 1,SYM, 1 .0 ) 

E0P=0.0 

DP 21 1=1,5 
IFII.EQ.5) FPP = 1 . 0 

CALI LRCURVIFRAMX, PRAM Y, 4,2, SYY,0. 0) 

21 CALL LRCURVI FR AMX I 5 ) , FR AMY ( 5 ) , 4, 2, S YM, EOP ) 
no 22 1=1,43 

22 CALL LRCURV(X,Y,0, l,SYM, 1.0) 

FOP=C.C 

DO 23 1=1,5 
IFII.EQ.5) FOP= 1.0 

CALL LRCURVI FRAMX , FRA YY, 4,2 , SYY.o .0 ) 

23 CALL LRCURVI FRAMX I 5) , F R AM Y ( 5 ) , 4 , 2 , S YM , E OP ) 
DO 24 1=1,5 

24 CALL LRCURVI X,Y,0, 1 ,SYM, 1.0 ) 

CALL L R M DF F 

RETURN 
ENTRY ot?NTI 

PRINTS INITIAL PARTICLE P'lS I T I n N S 


CALL LRMRG.N(XMRG,X«RG,YMRG, YMRG I 
CALL l RANGE I C . 0, XRNG.C.C , YRNG) 
CALL LRGRIDI 1,1, 0.0, 0.0) 

CALL LRXLFGICHARX.l ) 
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CALL LRY( EGICHARY, l) 

CALL LRCNVT (T , 3 , CHART ( 2 ) , 3, 5, 2 ) 

CALL LRTLFGI CHART, 12 > 

CALL LRCURV(BX,BY,5,2,SYM,Q.O) 

CALL L R NON 

CALL LRCURVI X,Y , MLAST ,5 , SYMBOL ,1 .0 ) 
CALL LRNOFF 
RETURN 
ENTRY PRNTN 

PRINTS NFVJ PARTICIE POSITIONS 


CALL I. R V R G N ( XMRG*XMRG*Y M RG»YMRG) 

CALL LRANGEI C. 0 , XRNG ,C . 0 , YR NG ) 

CALL LRGRir>< 1, l, 0.0, 0.0) 

YW = YW + DELT/DELY 

IF! YW.GE. FLOAT! II ) ■*- 0 . 4 9 ) YW=2.5 

CALL LRCNVT I T, 3 » CHART ( 2 ) , 3, 3, 2) 

CALL l.RCNVT ( SDMAX , 3,DE(3>,4,9,2) 

CALL LRLEGNIOE, 17,0, 3.0, 10.0,C.0) 

CALL LRCNVT ( CON MAX, 3,CF{ A) ,4,9,2) 

CALL LRLEGN ICE, 2 1,0, 5. 6,10.0,0.0 ) 

CALL LRCNVTI I TERO, 1, 10 I 3) , l ,4,0) 

CALL LRLFGNI 10, 12,0,3.5,0.4,0.0) 

CALL LRCNVTI ITER I, 1 , I K ( 3) ,1,4,0) 

CALL LRLFGNI IK, 12,0, 6.1,0. 4,0.0) 

CALL LRXLEGICHARX, 1 ) 

CALL LRYLTGICHARY, l) 

CALL LRTLFGICHART, 12 ) 

C POSITION OF MARKER INDICATING VELOCITY OF WALL 

CALL LRCURVI XW, YW, 1, 4, WL ,0. 0 ) 

C OUTLINT OF CAVITY 

CALL LRCURV(BX,BY,S,2,Sym,0.0> 

CALL LRNON 
CALL L R S ON 

CALL L RC UK V ( X , Y , ml AST , 5 , SYMBOL , l ,0 ) 

CALL LRSOFF 
CALL LRNOFF 
RE TURN 
ENTRY TAIL 
00 25 1 = 1,48 

25 CALL LRCURVI X ,Y ,0 , 1 ,SYM, 1 .o » 

RFTURN 
ENTRY PRNTF 

C PRINTS FINAL F R A m F OF MOVIE AN 0 R r ’S I T l '•«> S OF VO>TFX CENTER E 

C 

CALL LRMON 

CALL LRL EGNI TIN, 2’,C , 1 . 0 , 9 . 6 , I . ' ) 

C NORMALIZE THF VORTEX CENTER PCSfTIINS 

DO 2 1= 1 ,N1 AST 

VC XP ( I ) = W* ( 1 .0-1 VCX ( I l-2.5)/FJJ*?l 
2 VCYPI I ) = (VCY( I )-2.5)/FI I M2 
VCXN = VCXP ( ML AST ) 

Vf. YN = VC YP( NL AST ) 

CALI LRANGEIO. 0,H, 0 . 0 , W ) 

CALL LR OR 10(2,2, 0.5, 0.5) 

CALL LRMRGNI YMRG ,YMRG< XMRQ, XMRG ) 

CALL LRCHSZ ( 3) 

CALL LRTLFGI VXC, 63 ) 
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C POSITION OF UNSTEADY VORTEX CENTER AT END OF CALCULATION 

CALL LRCURVI VCYN.VCXN, l, 4, S VC U, O.O) 

C POSITION OF UNSTEADY VORTFX CENTER DURING CALCULATIONS 

CALL LRCHSZ(C) 

CALL l RCURVI VCYP, VCXP,NLAST,2,SYM,1.0) 

CALL LRGRID(2» 2 *0.25 *0.25) 

I VC = VC Y I NLAST) +0.5 
JVC = VCX I NL AST ) +0.5 
PL = 0.5 + VCY(NLAST) - FLOATIIVC) 

OL = 0.5 + VCXINLASTI - FLOAT(JVC) 

C TRAVFRSE PERPENDICULAR TO MOVING WALL 

DO A JJP=3,JJ 
JP = JJP -I 

VPROFLIJP) = PL*V( IVC, JJP,K ) + ( l . O-PL ) * V < I VC- l , JJP , K ) 

A VF RTPS ( J P) = W* ( 1.0- (FLOAT (JJP >-2 . 5 ) / FJ JM2 » 

CALL LRMRGM 1.0,0. A, 1.0,0. A! 

CALL LRANGEI-0.5, l.OtO.OfW) 

CALL LRXLFGI XVEL,4?> 

CALL LRYLEG ( YPOS, 17) 

CALL LRTLFG(TPAR,43) 

CALL l RCURVI VPROFL (2) , VE RTPS ( 2 ) , J JM? ,4, SVCU.0.0 I 
CALL LRCURV ( VPROFL, VER TPS, JJ,2, SC, 1.0) 

C TRAVERSE PARALLEL TO MOVING WALl 

DO 5 I l P = 3 , I I 
IP = IIP - 1 

UPROFL (IP) = -QL*U 1 1 IP , JVC » K ) - ( l .O-QL ) *U ( I I P , J VC- 1 , K ) 

5 HORZ PS ( I P ) = I FLOAT! I IP)-2. 5 )/FIIM2 
CALL LRANG6 (0.0, H,-0 .5, l .0 ) 

CALL LR XLEG ( XP0S.19) 

CALL LR Yl EG( YVEL,47) 

CALL LRT LEGIT PND, 50 ) 

CALL LRCURV ( HORZ PS ( 2 ) , UPROFL (2 ) , 1 1 M2 ,4 ,SVCU,0.0 ) 

CALL LRCURV ( HOR 7 PS, UPROFL , I I , ? , SC ♦ 1 . 0 ) 

C 3D PRESSURE PLOT 

IF(W.GT.H) JJM2=JJM2/2 
IF(W.LT.H) IIM2MIM 2/2 

CALL PL0T3D( 3.0 ,FJJ,3.0 , FI I , S , II M2, J J MR, [ l M ?, J JM 2, FPR F SS ,C UBE ) 
CALL ROTATE ( C. 0 ,0. 0 ,45. C ,. TRUE. , .FALSE. ) 

CALL ROT ATE1C.0, 35 .(' , G . 0 , . TRUF . , .F AL SF. ) 

RFTURN 

END 

CALL L R MON 
SUBROUTINE DANDR 

C CALCULATES DISCREPANCY AND SOURCE TERM FOR POISSON’S FQUATION 

COMMON D ELX, DE LY f OFLT,W,H, TLA ST, VISC, OMEGA, GX, GY, 

1 DEL XSQ , DEL Y SO , DX4 ,D Y4 , C. XQE L X ,GYDF L Y , V I SCX, VI SCY ,C OMOG , L1MGO 7 , 

2 XOY,XOYS, VX,VY , VXY, FDELX, FDEl Y, TDXDY, DXOD Y, DYOf) X, GVHGXW ,NL A ST , 

3 K ,L , 1 1 , JJ , I I PI , JJPl , I I Ml , JJM1 ,N,T .DELPHI , D M AX , CCNMAX, XM AX.NPART 
COMMON /U VP /U( 44,44,2) , V( 44,44, 2) ,PHI < 44,44) ,D( 44,44) ,R (44 ,44 ) , 

1 DCON » PCCN , M AX IT , I TERC , l T ER I , SQM AX 
DMA X= 0.0 
DO 22 J= 3, J J 
DO 22 1=3,11 

D( I , J ) = IlKI ,J,L»-U(I ,J-1,L) ) / DELX + (V( I , J,L)-V ( l-l , J,L ) ) / DELY 

IF (ABS(D( I,JI J.LE.OMAX) GO TO 22 
OMAX=ABS(D( I , J ) ) 

SDMAX=D( l , J) 
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C GENERAL TER M 

?2 R(I,J)= ((U( t,J+l,L)*U( I.JtL) )**2 - 2.0*(U( I ,J ,L )*U( I, J- l,L ))**2 
1 + ( UU » J— 1 1 U + IJ( I , J-2 .L) )**2 ) /DX4 

? * ( (V< I+l,J,L)*V( I,J,L>)**2 - 2.0*1 V< I ,J,L) +V t I-l , J ,L ) )**2 

3 + ( V< l-l . J,L) *V< I-?,J,L)J**?)/DY4 

4 + (UI(IM,J,LI»-UU,J,UI*(V!I,Jfl,LHV(!,J,L)l 

5 + (U< I , J-1,L) +U( l -1, J-l ,L H *< VI I-l, J,L ) + V( l-t , J-l ,L») 
ft ~ UJU , J,L)*IJ( I-l ,J,UI*(VI I-lt J + 1,L )+V< I-l.J »L ) ) 

7 - (U( I*1,J-1»L)+U(I »J— 1»H ) * ( V( I , J « L ) *V ( I tJ-l iL) ) ) / TDX DY 

8 - D(I , J ) / DEL T 
27 00 27 J = 3 * J J 

C ROW ABOVE UPPER WALL 

R ( 2 » J I = 0 . 25 * ( DYOOX*< ( U ( 3 , J , L ) +U { 2, J , L J) *( V< 2, J* l , L) ♦ V( 2 , J ,L ) ) 

1 - <U(3t J-l ,L> *U<2, J-l ,L) )*( V(2t J,L)+V(2, J-l, Llll 

2 + ( (V(3, J,L»*V(2,J,Ln**2-(V(2,J,L)*V( 1,J,L>)**2 )) 

3 - 0.5*VISCX*< r)YODX*(V(2,J+l,L)-2.0*V(?, J,L >+V< ?,J-l,L> > 

A- ( U< 3, J,L) -U< 3 t J-l,L) -U( 2 , J,L> + U(2 » J-l , L) >) -GYOELY 

C ROW BELOW LOWER WALL 

23 R ( 1 1 *1 * J I = -C.25 *(<( V( I I + lt J,L)+V( I I, J,L ) )**? 

1 - ( V( I l ,J,L)*V< II-1,J,L)I**2> 

2 + DYODX*( (U( I 1*1, J ,L ) +U( II , J,U ) *( V( 1 1, J + 1,L )* V( 1 1 , J ,L)) 

3 - (U(I 1+1 , J-l ,L)+U( l I , J-l ,L ) )*( V( l I , J ,L >+V< I l , J-1,L >>) ) 

4 * 0 « 5*V I SC X*( DYOQX*( VI II ,J+1 ,L>-2. 0*V( I l , J ,L) *V( II , J-l ,L > J 

5 - (U( I I +1 , J,L )-U( I I *1 , J-1,L >-U( I l , J,L ) +U< l I , J-i ,L )>) + GYOELY 
DO 24 1=3,11 

C COLUMN BEFORE LEFT WALL 

R ( I , 2 ) = 0.25*{ ( (U( I,3,L > +U ( I,2,LI)**2 - ( U ( I , 2 ,L > + U< l , 1 , L ) ) ** 2 ) 

1 + DXODY*( (ll( I +1 ,2 ,L) +U(I ,2 ,L ) )* ( V< I ,3 ,L) +V U ,2 ,L ) ) 

2 - (U( 1-1,2, Ll+Ul I * 2 * L ( )*(V(I— l,3,L)+V(I-l*2,L>) )) 

3 - 0.5*VISCY*< DXODY*< IJ(I *1,2,L)-2.0*U(I»?,L)*U( I - 1 , 2 , L ) ) 

4- ( V< I ,3,L>— V( I-l ,7,L)-V( I , 2 , L ) *V ( I-l, 2, LI) ) -GXDELX 

C COLUMN after right wall 

24 R ( I , J J* 1 ) = -0.25*1 I (U( I »JJ + 1 ,L )+U( I , J J , L I ) **2 

1 - (U( I , JJ,L! + U( I , JJ-l,l ))**2 ) 

2 + DXODY* ( (U< 1+1, JJ,L > *U( I, J J,L ) )*( V( I,J J+l ,L)*V< I , JJ,Lll 

3 - (U< I-l, JJ ,L ) *U( l , JJ,L ))*< V< l-l , JJ + 1 ,L) *V< I- 1, JJ,L > >)) 

4 + 0.5*VISCY*( OXODY*( U( 1* l , J J , L ) -2 . 0*U( I ,JJ,L» +UU-1 ,JJ,U) 

5 - (V< l, JJ+l.L ) -V ( I-L,JJ + 1,L )-V( I,JJ,L >*V( I-l, JJ,L>)> 

6 + GXDFLX 
RFTURN 

END 


SUBROUT INE PRESS 

C CONTROLS ITERATION OF POISSON*S EQUATION AND DETERMINES CONVERGENCE 

COMMON DELX,OELY,DELT,W,H,TLAST,VISC ,OMEGA , GX , GY , 

1 OELXSQ, DELYSQ, DX4, DY4, GXDELX, GYOELY, VI SC X , V I SC Y ,C OMOG , OMOO Z , 

2 XOY ,XOYS,VX,VY »VXY ,FDELX, FOELY, TDXDY, DXODY, DYODX, GYHGXW , NL AST, 

3 K , L , I I , J J , I I P l , J J P 1 , I I M 1 , J J M 1 , N , T , OE L PH I , DMAX , CON MAX , XMAX ,NPART 
C0MM0N/UVP/U(44,44,2 ) , V ( 44, 44 , 2 I , P H I ( 44 , 44 J , 0 ( 44 , 44) , R ( 44 , 44 > , 

1 r>CON,PCCN,MAXIT,TTERO, I TER I ♦ SDMAX 
DIMENSION CON ( 44,44) , TE ST< 44 , 44 ) ,C ONST ( 44 , 44 ) 

IL = 2 + ( 1 1—2 )/2 
DO 36 J = 2 » J JPl 
DO 36 1 = 2, I IP1 
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36 CONST ( I , J ) = 0.?5*( (U ( l * J , L ) +11 ( 1 1 J- L t L > )**2 
1 + ( V(I f J,L)+V< )**2)+GYHGXW 

IF{ ITERO.NE.l I GO TO 31 
IF ( ITERI.F0.5 ) GO TO 31 
r TER I = l TER I - 1 

31 DO 35 I OUT - l f MAX I T 

MAXIT IS MAXIMUM number of iterations 

ITERO = TOUT 
DO 3? I I N= 1 , ITERI 
CALL ITFR 

32 CONTINUE 

DO 33 J = 2, J J P 1 
DO 33 1=2,11 PI 

33 TF ST ( I , J ) = PHI(I,J) 

CALL ITER 

CONMA X = 0 • 0 

TEST FOR CONVERGENCE 

DO 34 J = 3 i J J 
DO 34 1=3,11 

CON ( I , J ) =ARS { PHK I« J)-TEST( I *J) )/(ABS( PHI(I,J) ) 

1 ♦ ABS ( TEST < I , J ) ) + CONST! It J) ) 

I F ( CONI I t J ) .GT.CONMAX ) CDNM AX =CON ( I , J ) 

34 CONT INIJE 

I F ( CONMAX.LE.PCON ) GO TO 37 

35 CONTINUE 
ITERO = MAXIT 
RFTURN 

37 TST = 0 . 25* ( PH I Cl L , JJ )+ PH I ( 1 1 J JP l ) + PH I( I L + 1 , JJ ) +PH I ( IL + 1 , J J P 1) ) 
CHG = 1.0-TST 

NORM AL I Z F TO PHI* 1.0 AT CENTER OF WALL OPPOSI TE MOVING WALL 

DO 38 J=2 , J JPI 
DO 38 1=2, IIP1 
PHI { It J) = PHI C It J ) ♦ CHG 

38 CONTINUE 
RETURN 
END 


SUBROUTINE ITER 

COMMON DELX,DELY f DELT, W,H, TLA^T, VI SC ,OMEGA,GX,GY . 

1 DELXSO, DELYSQ, DX4, DY4 , GXDEL X , GYDFL Y , V I SC X , V I SC Y ,C 0*00 , 0 M G O l f 

2 XOYfXOYS, VX , VY ,VXY ,FDELX, FDFl Y f TDXDY, DX 0 DY , D^ 0 Q v , r;v H C X W , N L AST, 

3 K,L, II, JJ, I IP It JJPlt riMlt JJM1 ,N,Tt DEL PH I , DK*AX,CCf' 'MX , X *< Ax , NP AR T 
COMMON/ UVP/U ( 44, 44 , 2 ) , V ( 44 ,44,2)» D HM 44 , 44 ) »D( 44,44) , R { 4 4 , 4 4 ) f 

l DC ON, PC ON, MAX l T , I TERO, I TFR I , S D M A X 

ROW ABOVE UPPER WALL 

DO 41 J=3 , J J 

41 PHI(2,J)= PHI! 3,J)+R(2,J) 

COLUMN BEFORE LEFT WALL 

DO 43 1 = 3, I I 

PH 1(1,2)* PH I ( It3)+RUt2) 

DO 42 J= 3, J J 
GENERAL TERM 

42 PHI (I y J) - COMOG*PH I ( I , J ) +OMGC ? * ( ( PH I f I t Jf* l ) «-PH I ( I , J - 1 ) ) / DFL X SO 
1 + (PHI ( 1 + 1, J) +PHI ( I- It J ) ) /OELYSQ + R l l , J) ) 



C COLUMN AFTER RIGHT WALL 

43 PH I ( I * JJ4-1 ) = PHI ( I«JJ)+R(I tJJ+ll 

C ROW BELOW LOWER WALL 

DO 44 J=3,JJ 

44 PHMII + I,J)= PHKII.JI + RUI + ltJ) 

RETURN 

END 


SUBROUTINE UANDV 

C CALCULATES AXIAL AND NORMAL VELOCITIES 

COMMON OELX, DELY,DELT, M, H,TLAST,V l SC , OMEG A , G X, GY , 

1 DELXSQ,DELYSQ,DX4,DY4 ,G XDELX , GYDELY , V I SCX , V I S CY ,COMOG, OMGO Z, 

2 XOY , XOYS, VX,VY, VXY.FDELX, FDEL Y , TD XD Y, DXOD Y, DYCDX , GYHGXW ,NLAST, 

3 K.L.I I t JJ »I I PI » JJPlt I I Ml, JJMl,N,T,OELPHI, DM AX , CONMA X , XM AX , NPAR T 
C0MM0N/UVP/UI44 ,44 , 2 ) , V I 44 , 44 , 2 ) , PH I ( 44 ,44 ) , D ( 44 , 44 » , R ( 44 , 44 ) * 

1 DCON,PCON,MAX IT, ITERO, I TER I , SDMAX 
DO 52 1=3,11 
DO 52 J=3» JJ 
C GENERAL TERMS 

U( I » J,K) = U< I, J,L)*DEIT* tUUU, J,L »+U( !,J-1,L » )**?. 

1 - ( U( l , J+l,L)+U(I ,J,L) )**2» /FDELX 

2 ♦ (<U( I,J,L)+U( I-l,J,L>)*(V< I-WJ4-1,L)+VI 1-1, J,L) ) 

3 - ( U< 1+1, J ,L ) +U( I , J,L ) >*( V( I , J + l ,L) +V ( I , J, l )) )/ FDFLV 

4 + GX+IPHI ( I, J )-PHI( I , J+m /DELX 

5 + VY*(U(in, J,L)-?.0*U( I,J,L)+U( I-1,J,L) » 

6 - VXY*( V( I , J+L ,L) -V( I -1 , J+L ,L)-V( I , J,U«-V (1-1 , J ,L )>) 

V(I,J,K)= V( I, J,L >+DELT*(< (U( H-l,J-l,U+U( I , J-l,LM 

1 *( VC It J,L >+V< I , J-l ,U )-(IJ( l + l, J,L >+U< I, J ,L ) ) 

2 *( V( l,J+l,L)*V(I ,J,L> ) ) ZFDELX 

3 + ( (VI I»J,L )+V( I- l, J , L I ) ** 2— I VI I+1,J,L)*-V( I ,J,L))**?) /FDELY 

4 + GY+(PHl( I ,J)— PHI ( 1+1, JU/OELY 

5 + VX*( V( I , J«-l,L>-2. 0*V( I , J ,L H-VI I , J-l , L) ) 

6 - VXYMUI 1 + 1, J,L )-U( I , J ,L I-U( t + 1, J- 1 ,L ) +U< I , J-l ,L)» ) 

52 CONTINUE 

DO 53 J = 3 , J J 
C LOWER NO SLIP ROW 

53 VIII , J , K ) = 0.0 

C RIGHT MO SLIP WALL 

DO 51 1 = 3, I I 

51 U( I ,JJ,K)= 0.0 

C CALCULATE BOUNDARY VALUES 

DO 54 J=3,JJM1 

C ROW ABOVE UPPER NO SLIP WALL 

U( 2, J,K) = -U( 3, J,K) 

V( 1, J,K)= -V ( 3, J , K ) 

C ROW BELOW LOWER NO SLIP WALL 

U( II + l, J ,K )=-U(I I , J,K) 

54 V(II+1,J,K>= — V ( I I — 1»J»K I 
V( 1 » J J » K ) = — V I 3 , J J , K ) 

V( 11*1* JJ.K )= -V( I I— 1, J J ,K ) 

DO 55 1 = 3, I I Ml 

C COt (JMN BEFORE LEFT MOVING NO SLIP WALL 

U(I»1»K)= — U(I,3,K) 

VI I , 2 , K ) = 2.0-VI 1,3, KI 
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C COLUMN AFTER RIGHT NO SLIP WALL 

U( I, JJ + 1 ,K )= -U( I, JJ-ltK ) 

55 V( I , JJ + 1 « K ) =— V ( I , JJ,K) 

IJI lit 1 * K ) =— U ( II * 3 * K ) 

IJ( II*JJ+1»K)= -!J( l I, JJ- 1,K ) 

C UPPER LFFT CORNER 

U( 2 * 2 * K ) = U( 3,2,K» 

V< 2,2, K) = V ( 2, 3,K ) 

C LOWER LFFT CORNER 

U( I I +1 , 2 *K )= U( II, 2, K) 

C V( I I ,2 ,K |= VI! I ,3,K) 

C UPPER RIGHT CORNER 

U ( 2, JJ,K ) = -U( 3, JJ ,K) 

V(2, JJ+1 , K ) - V ( 2 , J J , K ) 

C LOWER RIGHT CORNER 

IJ( I I+l»JJ,K) = — U ( I I , J J , K ) 

V( I I , JJ + 1 * K ) = V ( I l , J J ,K ) 

RETURN 

END 


SUBROUTINE VTXCTR 

C CALCULATES THE POSITION OF THE VORTEX CENTER 

COMMON DELX, DEL Y,DFLT ,W,H,T LAST, V ISC, OMEGA, GX, GY, 

1 DEL XSQ , DE L YSQ » D X4 , DY4 » G XDE LX , GYOF L Y » V I SCX , V I S CY » C OMOG, OMGOZ, 

2 XQY , XO Y S * VX , V Y , VX Y , FDEL X , FDEL Y , TD XD Y, DXOO V, OYOO X »GYHG X W »NL AST , 

3 K ,L , I I , JJ , I IP1 , JJPl , I I Ml, JJM] ,M,T, OEL PHI , DM AX , CON MA X , XM AX , NPAR T 
C OMMON/'J VP /U( 44,44,2) , V ( 44 ,44 , 2 ) , PH I ( 44 ,44 ) ,0(44,44) ,R(44,44), 

1 DCON, PC ON, MAX IT, I TERO , I TER I , SRM A X 
COMMON/ VC /VCX ( 1500 ) ,VCY (1500 ), JCH(50 ),NPPC 
FUNCT ( A , B ) = AB S ( A / ( A— 8 ) ) 

DO 03 1 = 3, I I M 1 
DO 01 J=3,JJM1 

IF ( V( I , J,K ) *V( I , J+1,K ) . LF.O. C) GO TO 92 
01 CONTINUE 

VCX(N) = 0.0 
VCY(N) = 0.0 
RETURN 

92 JCH( [ 1= J 

93 CONTINUE 

J C H ( 1 I ) = J C H ( IIM1) 

00 P4 1=12,1 I M 1 
IHOLD = I 

JO = JCH( 1-1 ) 

J1 = JC H ( I ) 

J 2 = J C H ( I +1 ) 

1 F ( U( I ,J1 ,K)*U( 1+1 ,J2,K).LE.C.T ) GO TO 95 

94 CONTINUE 
VCX(N) =0.0 
VC Y ( N ) =0 .0 
RETURN 

95 I F ( ABS (U( IHOLD + 1, J2,KM .GT. ABS< U( IHOLD, J t ,K) ) ) GO TO 96 

c UPPER HALF CFLL 

X 1= FLOAT( J1 )+FlJNCT( V ( I HOLD , J 1 , K ) , V ( I HOL D , J l +1 , K ) ) 

X 2= FLOAT ( J2 ) +FUNCT ( V ( I HOI D + 1 , J 2 , K ) , V( IHOLD+ 1 , J 2+ l , K ) ) 
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VI = FI OAK l HOLD 1+0.5 
GO rn 97 

C LOWER HALF CELL 

96 Xl = FLOAT! JO+FUNCT! V ( I HOL 0-1 , JO . K > , V ( I HOL D-l * JO + l , K ) ) 

X 2= FLOAT ( J L > * FUNC T ( V< I HOLD , J 1 , K ) , V (IHOLD , J 1 + 1 , K ) ) 

Y1 = FLOAT ! IHOLOI-O.5 
07 I F ( J2.NE.Jl ) 00 TO 103 

I F ( AflS < V! IH0LQ»J1+1»K) ) .GE • ABS < V ( I HOLD , J 1 , K ) ) ) GO TO 100 

C LEFT HALF CELL 

on 98 i = n,i imi 
im= i 

l F ( UU+l, Jl+l.K )*>M l, Jl+l,K).LE.0.0 ) GO TO 99 

98 CONTINUE 
VCXIN )=C.O 
VC Y ( N ) =0 .0 
RETURN 

99 Y3= FLOAT! IHOLO H-FUNCTI IK I HOI 0 , J 1 , K ) , U ( I HOLO+ l , J1 , K 1 1 

Y4 = FLOAT! I M I+FUNCT ( U ( l N, J 1 + l , K ) ,U < IM + 1, J 1+ l» K ) ) 

X 3 = FLOAT! J 1 ) + 0. 5 

GO TO 108 

C RIGHT HALF CELL 

IOC DO 101 1=11, I I M l 

IM= I 

IF! U( I , J l — 1 ,K) *U( 1+ 1, Jl-l ,K) . LE.0.0 ) GO TO 10? 

1C l CONTINUE 

VC X ! M ) = 0.0 
VCY(N)= 0.0 
RETURN 

102 Y3= FLOAT! I N) +FUNCT ( U U *» J l-l , K > , U ( IM «■! , J 1- l , K ) > 

Y 4= FinAT! IHOLDM-FUNCT! U( I HOLD , Jl , K ) , U ( I HOL 0+ l , Jl , K 1 I 

X 3 = FLOAT ( J 1 ) - 0.5 

GO TO 1 C 8 
C J2.NF.J 1 

103 DO 104 I = I HOLD, I IN l 
(Ml = I 

IF! IMI ,J1,K)*U( l+l, J1,K1. LE.P.O ) GO TO 109 

104 CONTINUE 

VCX(N) = 0.0 
VCY(N) = 3.0 

RFTURN 

109 DO 105 I=ll,TIMl 
IM2 = I 

IF! U! I , J?,K)*!J( I+1,J2,K I.LF.0.0 ) GO TO 106 

105 CONTINUE 

VCXIN) = 0.0 

V C Y ! N ) = 0.0 
RFTURN 

1 C 6 IF! J2.GT.J 1 ) GO TO 107 

Y3 = F| OAT! IM? ) +FUNCT ( U ( I M2 , J 2 , K ) , IM IM 2 + 1 , J ? , K ) ) 

Y4= FLOAT! IM1J + FUNCT! IK I Ml , J 1 ,K ) ,U ( I MH-1 , Jl , K ) ) 

X3 = FLO AT ! J 2 ) + 0.5 

GO TO 108 

107 Y 3 = FLOAT! I wi ) 4-FUNCT! U(I M 1»J1»K)»U(I M 1+1,JI , K ) ) 

Y4= FLOAT! IM2H-FUNCT! U ( I M ?»J?*K) ,11! I M 2 + 1 , J ? , K ) ) 

X3 = FLOAT(Jl) + 0.5 

loe RS = ! Y 3 -Y 1 + { Y4-Y3)* ! X1-X3) ) /! 1 . 0- ( Y4-Y3 ) *( X2-X1 ) ) 

VCY < N ) = Yl + RS 
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VCX(N) = XI + (X2-X1)*RS 

RETURN 

END 


SUBROUT IMF PART 

c GENERATES INITIAL PARTICLE DISTRIBUTION 

COMMON DELX,0ELY,DELT, W,H,TLAST, VISC ,OMEG A , GX , GY , 

1 DELXSQ, DELYSOt 0X4, DY4, GXDEL X , GYDEL Y , V I SC X , V I SC Y ,C OMOG , OMGOZ , 

2 XOY f XOYS, VX,VY, VXY,FOELX,FDELY,TDXDY,DXODY, DYOOX, GYHGXW,NL AST, 

3 K,L, Ilf JJ, IIPU JJPltllML, JJMl ,N,T, DELPHI , DM AX f CONMAX , XM AX * NPART 
C0MM0N/UVP/UI44, 44, 2 ) , V ( 44 , 44 , 2 ) , P H I ( 44, 44 ) * D ( 44 * 44) * R ( 44,44) , 

1 DCON,PCON,MAXIT,ITERO,I TER I ,SOMAX 
COMMON/M I /M, ML A ST, XI 1600 ) ,Y( 1600) 

COMMON/ VC/ VCX ( 1500 ) ,VCY ( 1500 ), JCH( 50 ),NPPC 
I I M2 T-2 I 1-2) 

M=0 

IFINPPC.EQ.4) GO TO 62 

C ONE PARTICLE PER CELL ( OE LX=DEL Y ) 

DO 61 J = 3, J J 
XJ = J 

DO 61 1=3,11 
Y I = I 
M = M + 1 
X( M ) = X J 
61 Y ( M ) = Y I 
ML AST=M 
RF TURN 

c FOUR PARTICLES PER CELL ( DE L X=DEL Y ) 

6? DO 64 J=3 , J J 

XJ=FLOAT( JJ-0.25 
DO 6 3 1=3, II 
Y I =F LOAT ( I ) -0 .25 
M=M+1 
X { M ) = X J 
Y( M)=Yl 

X( M+l IM2T)=XJ+0. 5 

Y(M+I I M2T ) = Y I 

M = M+1 

X( M) = XJ 

Y ( M ) =Y I +0 • 5 

X( MM lM2T)=XJ+0.5 

63 Y( M + I IM2T)-YI^0*5 
M=M+ 1 1 M2T 

64 CONTINUE 
Ml AST=M 
RETURN 
END 
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SURRHIT I NF MOVE 

C CALCULATES PARTICLE MOVEMENT 

COMMON DFL X, DELY.DELT, W,H, TLA ST, VI SC , OMEGA , G X , GY , 

1 nELXSQ,DFLYSQ,DX4 , DY4 , GX DEL X , C-Y DEL Y , V I SC X , V I SC Y , CCMOG , UMG0 7 , 

2 XQY, XOYS, VX,VY, VXY.FDELX ,FOEt Y ,TOXDY, DXOOY, DYCOX !, GYHGXW ,NLAST , 

3 K ,L . 1 1 , JJ, I IP1, JJPl, I [M l, JJMl,N,T,OELPHI ,OMAX,CCN M AX,XMAX ,NPART 
COMM0N/UVP/U144 ,44 ,2 I , V ( 44 , 44 , 2 ) , PH I ( 44, 44 ) , D( 44 , 4 4 ) , R I 4 4, 4 4 ) , 

1 DC ON, PC ON » MAXI T,I TERO.ITERI ,SDMAX 
COMMON/ MI/M, ml AST, X( 1600 I, Y( 1600) 

F(A,P,C,D)= (l.O-P)*(l .C-Q ) * A + P*I1.0-Q)*B + Q*fl.O-P)*C + P*Q*D 

00 75 M=1,MLAST 

C PARTICLE IS IN I,J CELL 

1 = Y ( M I + 0.5 
J = X ( M I + 0.5 

C FRACTIONAL PARTS OF PARTICLE POSITIONS 

FR AC Y= YIM) - FLOAT! I ) 

FR ACX= XIMJ - FLOAT ( J ) 

P- FRACX *■ 0.5 
I F ( FRACY.GE.0.0 I GO TO 71 

C UPPER HALF CELL 

0 = -FRACY 

IF { T.EQ.3 ) GO TO 76 

!JM = FI U< I , J- 1,K I ,U( I , J ,K) ,U( 1-1, J-I,K> ,U( 1-1 , J ,K) ) 

GO TO 72 

C LOWER HALF CFLL 

71 Q = 1.0 - FRACY 

I F ( I.EQ. I I ) GO TO 77 

UM = Ft Ut I + 1,J-1,K! ,U« 1*1 ,J,K> ,1111 ,J-1 , XI ,UII ,J ,K) ) 

72 Q = 0.5 - FRACY 

IF! FRACX. GE. 0.0 » GO TO 73 

C LEFT HALP CELL 

P = 1.0 + FRACX 

IF! J.EQ.3 ) GO TO 73 

VM = FI VII , J- 1,K) ,V( I , J ,K) , VI I-lt J-l ,K) , VI I-l , J ,K) I 
GO TO 74 

C RIGHT HALF CELL 

73 P = FRACX 

I F I J.EQ.JJ I GO TO 79 

VM = FI VII , J,K> ,VI I , J+l ,K> , VI I-l ,J ,K ) ,V I I-l , J+l ,K ) I 

74 X(M) = X I M ) + UM*DELT/DELX 

Y(M) = Y ( M » + VM*OEL T/DELY 

75 CONTINUE 
RETURN 

76 UM = (1 .0-P!*U(3, J-1,K ) + P*UI3,J,K) 

GO TO 72 

77 UM = I l . 0-P ) *U I IltJ-lfK) + P*ll( II * J , K> 

GO TO 72 

78 VM = ( 1.0-Q)*V? T ,3,K) + G* V I I -1 , 3 , K ) 

GO TO 74 

79 VM = 11 ,C-Q»*V( I, JJ,K> * Q*V ( I-l , JJ, K ) 

GO TO 74 

ENO 
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Bl PCK DATA 

-SFTS INITIAL CONDITIONS 

CO M MPN/UVP/U< 44, 44 ,?) ,V( 44 ,44 ,? 1, PH 1(44,44 ), 0(44 ,44) ,R (44,44 ) , 
l DCON, PC ON, WAX! T, I TERO, I TER I , SDMAX 
COMMCN/PL/VPRDFL (50 ) fUPROFL ( 5C ), VFRTPS ( 50 ) , HORZP S( 50 ) 

DATA U/ 3372*0.0/, V/45*0. C, 41* 2.0 ,1395*0.0,41*2 .0,1850*0.0/ 

OATA PHI/ L9 36* 1 . 0 / 

DATA DCON/3.5F-3/ ,PCON/? .0E-4/,MAXTT/100/ 

OA TA [ TFPD/100/, I TER! /I 5/ 

DATA V ER TPS / 50 *0 • C / , VPRO FL /5C* D • 0/ » HOR Z p S/5f*l.O/, - 
I UPR0FL/5C*0 .0/ 

END 
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Lewis Motion Picture (C-271) is available on loan. Requests will be filled in the 
order received. 

This tutorial motion picture "Computer- Generated Flow Visualization Motion 
Pictures" (12 min, color, sound) shows how a computer solution of the equations of 
fluid motion can be presented as a calculated flow visualization experiment. A motion 
picture of the flow is made by using special particles that are moved with the fluid, 
displayed on a cathode ray tube, and photographed with a motion picture camera. 

Lewis Motion Picture C-271 is available on request to: 

Chief, Management Services Division (5-5) 

National Aeronautics and Space Administration 
Lewis Research Center 
21000 Brookpark Road 
Cleveland, Ohio 44135 


CUT 

Date 

Please send, on loan, copy of Lewis Film C-271. 
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Street Number 


City and State 


Zip Code 


National Aeronautics and Space Administration 
Washington, D. C. 20546 
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07U 001 33 51 3CS 71.10 00903 
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If Undeliverable (Section 158 
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"The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the tvidest practicable and appropriate dissemination 
of information concerning its activities and the results thereof .” 

— National Aeronautics and Space Act of 1958 
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